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[i] We present a method with which we determined the local lunar gravity field model 
over the South Pole-Aitken (SPA) basin on the farside of the Moon by estimating 
adjustments to a global lunar gravity field model using SELENE tracking data. Our 
adjustments are expressed in localized functions concentrated over the SPA region in a 
spherical cap with a radius of 45° centered at (191. 1°E, 53.2°S), and the resolution is 
equivalent to a 150th degree and order spherical harmonics expansion. The new solution 
over SPA was used in several applications of geophysical analysis. It shows an increased 
correlation with high-resolution lunar topography in the frequency band / = 40-70, 
and admittance values are slightly different and more leveled when compared to other, 
global gravity field models using the same data. The adjustments expressed in free-air 
anomalies and differences in Bouguer anomalies between the local solution and the a priori 
global solution correlate with topographic surface features. The Moho structure beneath 
the SPA basin is slightly modified in our solution, most notably at the southern rim of the 
Apollo basin and around the Zeeman crater. 
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1. Introduction 

[ 2 ] The gravity field of a planetary body is often expressed 
in globally supported spherical harmonics, since these 
functions arise naturally from Laplace’s equation expressed 
in spherical coordinates. Their global support however 
means that a large expansion in spherical harmonics is 
needed to describe small-scale stmctures, which results in 
the need to estimate a large amount of coefficients from data. 
Depending on the data density which may vary geographi- 
cally, this is not always possible or stable. In those cases, 
either smoothed global solutions are obtained by applying 
regularization of the ill-posed inverse problem in the form of 
a Kaula mle [ Kaula , 1966] (which constrains each coeffi- 
cient to zero with a given uncertainty, thus suppressing 
spurious power in especially the higher degrees), or, alter- 
natively, different representations that describe the gravity 
signal locally, at smaller scales, can be used. A variety of 
representations exist and they have been applied success- 
fully to data from satellites orbiting many different planetary 
bodies. To list but a few examples: the use of mass con- 
centrations (mascons) to describe gravity on the Moon [Muller 


'RISE Project, National Astronomical Observatoiy of Japan, Oshu, Japan. 
2 CRESST/Planetary Geodynamics Laboratory, NASA Goddard Space 
Flight Center, Greenbelt, Maryland, USA. 

department of Physics, University of Maryland, Baltimore County, 
Baltimore, Maryland, USA. 

Copyright 2012 by the American Geophysical Union. 

0 148-0227/12/20 1 1 JE003 83 1 


and Sjogren , 1968], or on Earth using recent GRACE data 
[Rowlands et al ., 2010], analysis of line-of-sight accelera- 
tions or Doppler residuals for Venus using gravity distur- 
bances [e.g., Barriot and Balmino , 1992], and estimating mass 
anomalies on Jupiter’s moon Ganymede [Anderson et al , 
2004; Palguta et al , 2006]. Furthermore, McKenzie and 
Nimmo [1997] introduced a method to use line-of-sight 
accelerations to directly estimate the admittance between 
gravity and topography locally, from which elastic thickness 
estimates can be derived. This method was applied to lunar 
data by Crosby and McKenzie [2005]. 

[ 3 ] The lunar gravity field seems especially suitable for 
using local methods, because the density of available track- 
ing data has been extremely asymmetrical: standard 2-way 
and 3 -way tracking from Earth stations (where the uplink 
and downlink stations are either the same or different, 
respectively) leaves a gap in tracking over the farside of the 
Moon due to the 1 : 1 spin-orbit resonance of the Earth-Moon 
system. The Lunar Prospector (LP) mission (1998-1999 
l Binder , 1998]) flew at altitudes as low as 15 km in its 
extended mission, resulting in an excellent coverage of the 
nearside with high-quality tracking data. Despite the farside 
gap, high resolution models with a spherical harmonics 
expansion up to a maximum degree of 1 65 were determined 
from the tracking data by making use of the aforementioned 
Kaula rule [Konopliv et al , 2001]. However, it was also 
realized that the full signal in the tracking data could not be 
accounted for even with these large expansions, as signifi- 
cant signatures were reported to be seen in the data residuals 
[Konopliv et al , 2001]. This has lead to several recent 
analyses using local representations of the gravity field over 
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the nearside of the Moon: Sugano and Heki [2004a, 2004b] 
used line-of-sight accelerations to estimate mass anomalies 
on the nearside of the Moon, while Goossens et al. [2005a, 
2005b] used Doppler data residuals directly to estimate 
gravity anomalies on the lunar surface. Again using line-of- 
sight accelerations, Han [2008] employed localized basis 
functions to estimate regional gravity fields on the nearside 
of the Moon for four locations, up to a resolution in spherical 
harmonics comparable to degree and order 200, and recently 
Han et al [201 1] updated this analysis to include the whole 
of the nearside at the same resolution. 

[ 4 ] The farside gap in tracking data was only recently 
largely filled by the Japanese SELENE mission [Kato et al , 
2010], which employed 4-way Doppler data [Namiki et al, 
1999] between a relay satellite in a highly elliptical (polar) 
orbit and a main orbiter in a circular (also polar) orbit at an 
average altitude of 100 km above lunar surface while the 
latter was over the farside. This has resulted in lunar gravity 
field models expressed in spherical harmonics up to a max- 
imum degree and order of 100 [Namiki et al, 2009; 
Matsumoto et al, 2010]. These models improved the view of 
farside gravity dramatically, revealing farside free-air grav- 
ity anomalies of large impact basins as negative rings with 
central highs, unlike their counterparts on the nearside of the 
Moon. 

[ 5 ] In this paper, we have used the farside Doppler data 
from SELENE to investigate whether we can extract more 
information from them than can be modeled with global 
spherical harmonics. The apolune of Rstar, the relay satellite 
involved in the 4-way tracking, was located over the south- 
ern hemisphere of the farside, resulting in a denser data 
coverage there [Matsumoto et al, 2010]. Figure 1 shows the 
4-way residuals with respect to the SELENE lunar gravity 
field model SGMIOOh [Matsumoto et al, 2010], together 
with the topography as measured with the laser altimeter 
(LALT) on-board SELENE [Araki et al, 2009]. A model of 
degree and order 100 in spherical harmonics already pushes 
the limits of what can be resolved from an altitude of 100 km 
above lunar surface, but the data coverage and remaining 
residuals as shown in Figure 1 indicate that there is still 
something to be gained from a localized analysis, either in 
terms of obtaining an adjustment with higher resolution, or 
in terms of redistributing the information in the spectral 
domain. In addition, the southern hemisphere of the Moon is 
dominated by the South Pole-Aitken (SPA) basin, the largest 
impact basin in the solar system (or second largest if the 
northern lowlands of Mars are counted as an impact basin). 
A localized analysis of this area might provide information 
for studying the structure of large impact basins, especially 
at smaller scales inside the basin. 

[6] For our analysis, we express the gravitational potential 
with localized spherical harmonic functions as described by 
Han [2008], which in turn was based on localized analysis 
on the sphere as given by Wieczorek and Simons [2005], 
Simons et al [2006], and Simons and Dahlen [2006]. Using 
SELENE tracking data over only the area of interest, we 
estimate adjustments to the coefficients of these localized 
functions, rather than estimating adjustments to the coeffi- 
cients of a spherical harmonic expansion directly. We com- 
pare our localized solution to solutions obtained from 
common global spherical harmonics, and we evaluate the 


performance of our solution in geophysical analysis using 
high resolution topography of the Moon. 

[ 7 ] This paper is structured as follows. We first introduce 
the data and methods in section 2, and in section 3 we 
include a benchmark test with LP extended mission data to 
establish the validity of the method and its implementation. 
In section 4 we show the results we obtained for the SPA 
area. The performance of the method is discussed in section 5, 
followed by the conclusions in section 6. 

2. Data and Methods 

[8] The tracking data we use will be the data from 
SELENE (apart from our benchmark test case using LP data 
as described in section 3). The main interest is in the SPA 
area, which means that the focus will be on using SELENE 
4-way data. However, 2-way data are also available, because 
a satellite can still be tracked from Earth beyond the limbs 
and poles of the Moon because of lunar librations. For 
instance, Figure 1 showed 4-way residuals with several gaps 
close to the south pole, which will be filled with 2-way data. 
This also means that there is a lot of LP data over the south 
pole, but for reasons explained in section 5 we did not 
include them in our local analysis. 

[ 9 ] The SELENE Doppler tracking data have a count 
interval of 10 s. In order to try and obtain more high- 
resolution information from them, it is possible to use shorter 
count intervals, but since that will also increase the noise 
level on the data, we did not pursue this here. Tracking data 
are processed with the GEODYN II software [Pavlis et al, 
2006]. While acceleration data as derived from Doppler 
residuals have proven to be a very powerful data type for 
local analysis, we chose for an integrated approach, where 
the local adjustments are estimated from the tracking data 
residuals directly. This choice was driven by the fact that the 
SELENE 4-way data are integrated data, i.e., they consist of 
the accumulated Doppler shift between the various links. 
There was no time-tagging for the separate links from which 
one could isolate only the Doppler shift between relay and 
main satellite. Our processing of the tracking data, in terms 
of precision force and measurement modeling, as well as that 
in terms of which parameters are estimated, is exactly the 
same as that described by Matsumoto et al [2010]. The 
difference between the approach described there and the one 
adopted in this work lies in the way that the gravity para- 
meters are estimated. 

[ 10 ] Our localized representation of the gravity potential is 
based on the work of Han [2008] (and references therein), 
where Slepian functions were used to estimate the gravity 
field over certain areas on the nearside of the Moon from 
LP line-of-sight accelerations. The Slepian functions are 
formed from linear combinations of spherical harmonics. 
For clarity, we repeat the basic transformations between the 
two sets here, and refer to the aforementioned literature for 
the details. 

[ 11 ] A function f{6, A) on the sphere is expressed in 
spherical harmonics as 

L l 

/(m) = Y E U-lm Y i m ( 1 ) 

1=0 m=—l 
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Figure 1. Residuals of 4-way tracking data with respect to the SGMIOOh model, with the focus on the 
SPA region (indicated by a red circle with a spherical radius of 45° centered at 191. 1°E, 53.2°S). Topog- 
raphy from SELENE altimetry is also shown. 


where 6 is the co-latitude, A the longitude, / and m are the 
degree and order, respectively, L is the maximum degree of 
the expansion, ui m are the coefficients, and F/ m (0, A) are 
the normalized spherical harmonics functions. This same 
function f{6 , A) can be expressed in the alternative basis of 
Slepian functions as 


for this expansion can then be found from the Slepian 
coefficients g n s /m and the spherical harmonics coefficients 
Ui m following \Han, 2008] 

l i 

V n = ^2 ^2 Uimgnjm (4) 

/=0 m=-l 


(L+ 1) 2 

AO, A)= Y V ngn{e,\) 

n = 1 


( 2 ) 


with different coefficients v n and Slepian basis functions 
g n {6 , A) [e.g., Simons and Dahlen, 2006; Han , 2008]. These 
new basis functions g n (6. A) can also be expressed in a 
spherical harmonics expansion 

L l 

g n (ff,A) = £; Y,snMe,\)- (3) 

/=0 m=—l 

[ 12 ] Equation (3) shows the transformation between the 
bases of spherical harmonics and Slepian functions. To 
achieve the desired localization of the Slepian functions, 
the coefficients g„ 5 /m are determined from maximizing the 
ratio between the energy within the desired area, and the 
energy over the entire sphere (this ratio is also denoted as 
concentration factor). As discussed by Wieczorek and 
Simons [2005], this can be posed as an eigenvalue prob- 
lem, where the concentration factors are obtained from the 
resulting eigenvalues, and the Slepian coefficients g„ 5 /m are 
given by the accompanying eigenvectors. Finally, in the 
case of the function^#, A) of equation (2), the coefficients v n 


This results in two equivalent expressions in different 
basis functions. 

[ 13 ] Han [2008] then separates the spherical harmonic 
functions into two parts uf m and u^n (or, equivalently, u\ and 
u 2 , as we will call them later), with the former consisting of 
those spherical harmonics concentrated in the desired area 
up to a chosen concentration factor y 0 (i.e., the concentration 
factor for the ftth Slepian function satisfies > 70 ), and the 
latter consisting of the remaining functions outside the cho- 
sen area. The full gravitational potential uj m satisfies 
Uim = uf m + Ufa. Following equation (4), the spherical har- 
monics expansion can thus be transformed into Slepian 
functions using a linear transformation u = Gv, with u 
consisting of the vector of spherical harmonic coefficients 
uj m , v consisting of the vector of Slepian coefficients v n , and 
G the transformation matrix between them. By splitting this 
into the two areas according to: 


u — ui + U2 — [Gi G2] 


vi 

V 2 


(5) 


Han [2008] then estimated only the Slepian coefficients Vi, 
i.e., only those that have their concentration in the desired 
area, using line-of-sight accelerations as the data type. 
This work was further expanded to global lunar gravity 
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field modeling from dynamically created normal equations 
(see below) by applying localized constraints based on this 
principle by Han et al. [2009], and these localized con- 
straints were applied to the full set of NASA Moon orbiter 
tracking data by Mazarico et al [2010]. In addition, Han 
et al [2011] recently used the method of Han [2008] to 
estimate the gravity field for the whole nearside of the 
Moon, again from line-of-sight accelerations. In our analy- 
sis, we follow the same principle of dividing the spherical 
harmonics coefficients into two parts: the main part that 
we are interested in is concentrated over the SPA area, and 
the remaining coefficients describe the potential outside 
of this area. 

[ 14 ] Because the focus is on a relatively small area on the 
Moon, long-wavelength signals might not be recovered or 
represented well in this method. This was noted by Han 
[2008], where it was shown that including the lower spher- 
ical harmonics degrees into the localized representation 
introduced an error of several tens of mGal (1 mGal = 
10 _5 m/s 2 ), whereas this error was much smaller when 
components with degrees higher than 30 were used, in the 
case of a spherical cap with a radius of 20°. This is thus an 
error of the representation form. We found similar issues in 
our local solutions, which are created from the localized a 
priori model plus the estimated adjustments. If the lower 
degrees of the a priori model were included in the localized a 
priori part of the constructed solution, we found a strong 
north-south component in the recovered anomalies that dis- 
appeared when the lower degrees were not included. We 
thus do not include the lower degrees of the a priori model in 
the localization. Instead, when we reconstruct the full gravity 
signal in standard spherical harmonics, we add the (global) a 
priori model up to / = 30 after the estimation. There is 
however one caveat when applying this remove-restore 
technique: the local basis functions have power over all the 
degrees, so the adjustment, when transformed back to stan- 
dard spherical harmonics, will show non-zero signal for 
/ < 30. We found however that the amplitude of the adjust- 
ments below / = 30 is small, at an RMS (root-mean square) 
of around 2 mGal, and the performance of the model as 
discussed later in sections 4.3 and 4.4 was unaffected. We 
also tried to mitigate the effect on the lower degrees in other 
ways, by generating the design matrix A with only entries 
for coefficients with / > 30, or by band-limiting the trans- 
formation matrix G, to span only the range 30 <1<L, which 
makes sure that coefficients with / < 30 do not get adjusted. 
The latter requires special attention: as pointed out by 
Wieczorek and Simons [2005] and Simons et al [2006], the 
eigenvectors needed to create the localized basis functions 
are most easily found from a tridiagonal matrix that com- 
mutes with the localization matrix, following results by 
Griinbaum et al [1982]. This works for the hill spherical 
harmonics range, but for the limited range 30 <1<L, one has 
to resort to a direct computation of the localization matrix, 
and determine the eigenvectors from this matrix rather than 
from the tridiagonal one. Wieczorek and Simons [2005] 
point out that this might be less stable, but we encountered 
no discernible differences. Both methods however (limiting 
the range of the partials, or that of the transformation matrix) 
did not lead to large differences between solutions. They 
were again at an RMS of a few mGal, and this shows that the 
residuals themselves as used in the adjustment already 


mostly have their signal confined to the high-frequency part 
of the gravity spectrum. We thus assume that the lower 
degrees of the gravity field model are already determined 
from a global rather than a local analysis. 

[ 15 ] We use a dynamical approach to process the tracking 
data. The data are divided into continuous time spans called 
arcs, during which the orbit is integrated using precision 
force modeling, and parameters valid for that arc only (such 
as initial state vector adjustments, measurement biases, non- 
conservative force model coefficients) are adjusted. Once 
these are converged (typically when the change in the 
weighted root-mean square of data residuals, consisting of 
the difference between the observations and those computed 
from precision modeling, is less than 2% between itera- 
tions), we generate the design matrix for the data used in the 
arc. The design matrix describes how the data depend on the 
model coefficients, and in standard linear form this is written 
as y = A u + e. Here, the vector y contains the data residuals, 
u is the vector of adjustments to the aforementioned spher- 
ical harmonics coefficients that are to be estimated, e are the 
data errors, and A is the design matrix that contains the 
partials of the observations with respect to the estimation 
parameters. This estimation problem is often solved using 
least-squares, which results in the need to form normal 
matrices, which are given as A r C^ ! A, where T denotes the 
transposed of a matrix, -1 the inverse, and C d is an a priori 
covariance matrix for the data. 

[ 16 ] For standard global gravity field determination, this 
normal matrix is formed from the complete design matrix A 
containing all the data included in the arc. In our localized 
approach however, we use all the data to converge the arcs, 
and then we select only the data over the area of interest for 
the formation of a reduced design matrix A red . We use only 
the data over the area of interest because in general those 
residuals are mostly affected by the gravity of features 
directly beneath them (depending on the data type used). 
This does mean that at the boundary the solution will be 
affected by features just outside of the area of interest. We 
discuss this in sections 3 and 4.3. We then use this reduced 
design matrix to form a normal matrix, and then proceed 
as usual by aggregating the normal matrices of all arcs 
involved, using the SOLVE software [ Ullman , 2002]. In the 
case that there are both arc and common parameters (those 
valid and constant throughout all arcs, such as the gravity 
field coefficients) present, the arc parameters can be elimi- 
nated through partitioning [e.g., Kaula , 1966]. We then 
transform the combined normal matrix following equation 
(5), and estimate only the coefficients of \ x which are con- 
centrated over the same area as where we collected the data. 
This is akin to a truncated singular value decomposition 
approach [e.g., Menke , 1989], where only those coefficients 
that are observable from the data (i.e., having non-zero sin- 
gular values) are estimated. 

[ 17 ] As pointed out by Simons et al [2006], and stressed 
also by Han [2008], the use of Slepian functions that are 
concentrated over a certain area requires far less coefficients 
than when a full spherical harmonics expansion is used to 
obtain the same resolution. This also means that less coef- 
ficients need to be adjusted from the data, which in general 
makes the inverse problem more stable, provided of course 
that the data actually have the information that is to be 
extracted. By choosing only data over the area, we also 
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reduce the number of computations that are needed to form 
the normal matrices per arc. In this way one is able to 
compute local solutions with high resolution much faster 
than global solutions with a similar resolution. Another 
advantage of the Slepian functions is that they can be 
transformed back and forth between spherical harmonics, 
which makes them easy to use with existing tracking data 
processing tools that are based on a spherical harmonics 
expansion of the gravity field. 

3. Benchmark Test With LP Data Over 
Mare Serenitatis 

[is] In order to test the suggested method and its imple- 
mentation, we created a benchmark test for a case where the 
gravity field of the Moon is well-known: over the nearside 
of the Moon. We chose to test the method over the Mare 
Serenitatis area, using the same approach as was done by 
Goossens et al [2005b], where the Lunar Prospector 
extended mission data over this area were used for local 
gravity adjustments. LP extended mission data were col- 
lected from January until July 1 999, at an average spacecraft 
altitude of 30 km [Konopliv et al., 2001]. These data thus 
contain a lot of gravity field information at the short wave- 
lengths. By using a pre-LP gravity field model in the pro- 
cessing, we can compare our local solution with known 
high-resolution global LP models, and thus we can infer 
whether or not the suggested method can extract the infor- 
mation that is in the tracking data. Because we use real data 
directly rather than conducting simulations, we can not 
quantify the recovery error, but this does give us the 
advantage of directly knowing how the method performs in 
a real-data environment, while the use of the older gravity 
field model ensures that we can also test whether the gravity 
field can be recovered at high resolution. We also do a 
similar test for our area of interest with the SELENE data, 
described in section 5.2, although that test is more meant to 
probe the influence of the a priori model used, rather than to 
serve as a benchmark test. 

[ 19 ] We used the same short-arc approach as used by 
Goossens et al. [2005b], where LP extended mission track- 
ing data were processed in short batches, when the satellite 
was over the Serenitatis area. This leads to arcs of typically 
ten minutes. We start with an estimate of the initial state 
that has been determined from longer arcs, using a high- 
resolution LP gravity model (LP150Q, a follow-up to 
the 165th degree and order model called LP165P described 
by Konopliv et al [2001]; the LP150Q model can be 
found through http://pds-geosciences.wustl.edu/missions/ 
lunarp/shadr.html). Starting from this estimate minimizes 
initial orbit errors. We then process the short-arcs with the 
pre-LP GLGM-2 lunar gravity field model [ Lemoine et al , 
1997] as a priori model. This model shows a much 
smoother signal for Mare Serenitatis than the LP models. 
While processing the tracking data with GLGM-2, we fixed 
some of the Keplerian orbit elements in order to ensure 
convergence of the short-arc orbit determination problem. 
Fixed parameters were the same as those of Goossens et al 
[2005b]: the eccentricity, ascending node, and argument 
of pericenter. We also noted that the range data are not 
as sensitive to local gravity variations as the Doppler 
data, since including the range data distorted the solutions. 


Doppler data have a better sensitivity, while derived (line-of- 
sight) accelerations as applied in many examples (mentioned 
in the introduction) are likely to be more superior, as is 
especially shown in the recent results of Han et al [2011], 
where high correlations between gravity and topography are 
achieved from line-of-sight acceleration data. 

[ 20 ] We estimate gravity field coefficients equivalent to 
a spherical harmonics resolution of degree and order 100, 
using Slepian functions concentrated in a spherical cap 
of 20° radius, centered around (18°E, 25°N). As described 
above in section 2, we do not adjust coefficients up to degree 
/ = 30. The used concentration factor was 0.0001, and this 
leads to having to estimate 477 Slepian coefficients. For 
comparison, a global spherical harmonics model up to degree 
and order 100 contains 10,197 coefficients. The gravity 
anomaly solution of Goossens et al [2005b] contained 806 
discrete anomalies in 1° x 1° blocks, over a square area 
(ranging from 5°-30° in longitude, and from 10°-40° in 
latitude). Tracking data residuals inside this area were used 
in the estimation process. Because the LP extended mission 
data were collected at a low altitude, downward continuation 
has only a small effect on the errors in the solution, and 
the local solutions can be obtained without any kind of reg- 
ularization to stabilize the inverse problem [Sugano and 
Heki , 2004a; Goossens et al., 2005a, 2005b; Han , 2008]. 

[ 21 ] Figure 2 shows gravity anomalies over the Serenitatis 
area, computed from either GLGM-2, LP150Q, or from the 
newly estimated model, with all expansions evaluated up to 
their respective maximum degree and order. It is immedi- 
ately clear that the local solution resembles that of LP150Q 
much more than it resembles the smooth GLGM-2 solution, 
as the increase in resolution inside the Serenitatis mascon is 
readily confirmed. The local solution also resolves the area 
around the mascon better than GLGM-2 does, when com- 
pared with LP150Q: the negative anomalies surrounding the 
mascon come out clearer, and they correlate well with the 
signal as seen in the anomalies from LP150Q. 

[ 22 ] Figure 3 shows the anomaly differences between the 
local solution and anomalies from either GLGM-2 or 
LP150Q. In addition, to assess boundary effects it also 
shows the difference between two local solutions with the 
functions concentrated in a spherical cap of either 20° or 
10°. The differences with anomalies from GLGM-2 (which 
are in fact the adjustments estimated from the data) espe- 
cially indicate that the extra signal from the data was 
extracted as there is clear structure visible in the differences. 
Differences with LP150Q show remaining signal as well, 
extending into the mascon itself, although most differences 
are concentrated at the boundaries. This is mostly due to 
boundary effects that are often present in local solutions, and 
which arise because although data are mostly affected by the 
gravity of features directly beneath them (as stated earlier), 
they also carry information from the surrounding area. 
Because we used the same processing (including the same 
way to generate the tracking data residuals) as Goossens 
et al [2005b], which used spacecraft tracks covering the 
square area listed above, the data do not fully cover the 
chosen spherical cap area and this exacerbates the boundary 
effects. In order to assess the influence of data coverage, 
we also created a solution that estimated only coefficients 
in a spherical cap of 10°. This solution is still sensitive 
to boundary effects (data outside the cap are also sensitive 
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Figure 2. Free-air gravity anomalies over the Serenitatis area computed from GLGM-2, LP 1 50Q, and the 
local solution. Each solution is expanded up to its maximum resolution in spherical harmonics, thus 70, 
150, and 100, respectively. 


to gravity features outside the area of interest that are not 
accounted for), but a comparison between the two solutions 
indicates where they are close to each other, and thus where 
boundary effects are smaller. The difference plot included in 
Figure 3 shows that in general the differences inside the 10° 
area are much smaller than the adjustments themselves, but 
they still show signal in what appears to be close to a zonal 
trend (positive and negative adjustments in the north-south 
direction). High-frequency differences with LP150Q can 
also be seen inside the mascon. These might stem from 
aliasing of higher degrees (note that the local solution pre- 
sented here is up to degree and order 100 while the LP data 
have information at higher resolutions), but on the other 
hand, these were also somewhat visible in the results of 


Goossens et al. [2005b]. Finally, correlations between vari- 
ous solutions over the mascon area were also computed, and 
it was found that the local solution is closer to LP150Q’s 
(the correlation between them is 94%) than GLGM-2’ s 
(a correlation with the local solution of 90%). The correla- 
tion between LP150Q and GLGM-2 is 89%. 

4. Results for SPA Using SELENE Tracking Data 

[ 23 ] The straightforward benchmark test served as a proof 
of concept, and it showed that the method as implemented, 
using Slepian functions, can indeed extract the extra infor- 
mation present in the tracking data. We now apply this same 
method to SELENE tracking data. 



q mGal 


- 125-100 -75 -50 -25 0 25 50 75 100 125 


Figure 3. Differences in free-air gravity anomalies between the local solution and either GLGM-2 or 
LP150Q. To assess boundary effects, the difference between two local solutions with the functions con- 
centrated in a spherical cap of either 20° or 1 0° are also plotted. 
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4.1. Processing of the SELENE Tracking Data 

[24] We concentrate on the South Pole-Aitken (SPA) 
basin. As mentioned in section 2, we process the tracking 
data in the same way as was done for the SGMIOOh model 
[Matsumoto et al., 2010], with the difference the a priori 
gravity field model that is being used (which was SGMIOOh 
itself in this analysis). Design matrices per arc are generated, 
including arc parameters such as initial state vector adjust- 
ments and measurement biases, and from them we select 
only data over the area of interest. All the parameters that 
are adjusted per arc are also included in the reduced design 
matrix. The arc lengths used varied per satellite [see also 
Matsumoto et al , 2010]: for Rstar they were on average 
2.6 days, while for the main satellite they were 12 hours (and 
6 hours after a reaction wheel failure in July 2008). The 
orbits for Rstar and the main satellite were determined 
simultaneously. The arc parameters are eliminated by parti- 
tioning when the normal matrices are accumulated. The 
influence of these parameters is discussed further in 
section 5.1. We used the following data weights in our 
solution: 2 mm/s for main orbiter 2-way Doppler data, 
1 mm/s for relay satellite Doppler data (due to different 
ground stations being used), and 1 mm/s for 4-way data. 

[25] We adopt the coordinates of the center of the SPA 
basin as given by Garrick-Bethell and Zuber [2009], i.e., 
(191. 1°E, 53.2°S). We include data in a spherical cap area 
with a radius of 45° around this center (this area is also 
indicated in Figure 1, which showed the SELENE 4-way 
residuals), and we also estimate the Slepian functions coef- 
ficients concentrated in this area, with a concentration factor 
of 10 -8 , and an equivalent spherical harmonics resolution of 
150 degrees and orders, which results in estimating 4,446 
Slepian function coefficients (compared to 22,797 coeffi- 
cients for a full spherical harmonics solution of the same 
resolution). This area is slightly larger than the SPA area 
itself, to account for boundary effects, which we discuss in 
section 4.3. 

[26] Considering that the tracking data are obtained from 
SELENE’s main orbiter in an orbit at 100 km above the 
lunar surface, this resolution is likely beyond the sensitivity 
of the data. We use this resolution because the use of Slepian 
functions allows us to estimate many fewer coefficients than 
when we would estimate an equivalent spherical harmonics 
model. Nevertheless, we found that we needed regulariza- 
tion for this data set to damp the amplitude in the estimated 
coefficients, despite the fact that local methods are often 
applied to circumvent the need for regularization. This need 
was to be expected however, mostly because of the main 
orbiter altitude at which the downward continuation of errors 
is substantial (as opposed to the LP extended mission results 
presented in the previous section, which were at a much 
lower altitude), but also from the amount of 4-way data 
available. For SELENE, 2-way data outnumbered 4-way 
data by a great amount [Matsumoto et al , 2010]. The total 
number of data into our solution was 203,927 points, of 
which roughly 170,000 were 2-way data, and only 34,000 
were 4-way data. The concentration factor (chosen to be 
10 -8 from an analysis of the representation error) affects the 
number of coefficients that will be estimated, but we found 
that changing this to 0.0001 (resulting in 3093 coefficients) 
did not take away the need to regularize the solution. 


Differences between solutions with these two concentration 
factors were small and located at the boundary. 

[27] To regularize the inverse problem, we constrained 
each coefficient toward zero with an a priori uncertainty 
following Kaula’s rule [ Kaula , 1966], 07 = /3// 2 , with / the 
spherical harmonics degree, 07 the uncertainty for the set of 
coefficients from the same degree, and /3 a constant. Since 
this is expressed in spherical harmonics, regularization was 
applied to the spherical harmonics normal matrix, after 
which the transformation following equation (5) was applied 
[see also Han et al , 2009]. The constant we used was 
3.6 • 10~ 4 , the same as was used for the SGMIOOh model 
[Matsumoto et al , 2010]. We also tested other constants and 
a power law based on the gravity potential from uncom- 
pensated topography, as was also done by Mazarico et al 
[2010], but we found that this did not influence the final 
results much, as the influence of the different regulariza- 
tion rules was confined to the higher degrees, where corre- 
lations with topography have dropped, see also results in 
section 4.3. 

4.2. Local Gravity Results for SPA 

[28] With the aforementioned processing, we obtained 
gravity adjustments over the SPA area as depicted in 
Figure 4 as free-air anomalies. The adjustments are in gen- 
eral small, between plus and minus 68 mGal, with an RMS 
of 17 mGal. Foremost, this indicates that the a priori model 
SGMIOOh already captures the information in the 4-way 
data well. The adjustments do show clear signals, with 
an outstanding circular anomaly signal which is located 
over the southern rim of the Apollo basin (location 36.1°S, 
151.8°W, diameter 537 km). There is also a distinct north- 
south pattern visible in the adjustments. This trend follows 
the ground track of the satellites involved, which is also 
seen in other high resolution lunar gravity field models, 
especially over the farside when there are no direct tracking 
data available [e.g., Konopliv et al , 2001]. In our case, 
we think this is mainly caused by the 2-way data. The 
line-of-sight geometry around the south pole is such that 
the 2-way data are likely to be less sensitive to radial gravity 
signals, since their line-of-sight is close to perpendicular 
to the radial direction there. For the 4-way data, the sensi- 
tivity depends on the exact geometry between the main 
orbiter and the relay satellite. The 4-way residuals are 
thought to describe mostly the remaining signal between 
these two, despite the data being integrated data over the 
various links. Because the relay satellite is in a higher 
orbit with the apolune over the lunar southern hemisphere, 
it should be less affected by gravity perturbations. This is 
also corroborated by small 2-way residuals for the relay 
satellite [see Matsumoto et al , 2010]. In Figure 5 we plotted 
the angle between the radial direction, defined as the direc- 
tion from the center of the Moon to the ground track point 
of the main orbiter in lunar-fixed coordinates, and the vector 
from the main satellite to the relay satellite. If this angle is 
small, the 4-way link between relay and main satellite is in 
line with the radial direction. Figure 5, however, shows that 
there are lots of tracks with a large angle, and they tend to 
correlate visually with the tracks in Figure 4. The histogram 
included in Figure 5 that shows the number of data points in 
angle bins of 5° also shows there are many data points with 


7 of 16 



E02005 


GOOSSENS ET AL.: SOUTH POLE-AITKEN LOCAL GRAVITY 


E02005 


Local adjustment Local solution 



I ~ " — ~ ~ = mGal I I mGal 

^ ' 1 T T T 1 ' 1 T ^ 1 T 1 T 1 ' 

-60 -40 -20 0 20 40 60 -400 -200 0 200 400 

Figure 4. Adjustments in (free-air) gravity anomalies for the SPA region from the local method, together 
with the full solution (a priori SGMIOOh and local adjustment). The projection is gnomonic, centered 
around 191. 1°E, 53.2°S. The names of several features on the lunar surface are also included in the adjust- 
ment plot, for reference. They are left-justified at the center of the feature, which means that the feature 
that is referred to is located at the start of its name. 


large angles. This means that some of the tracks as seen in 
Figure 4 might also be due to 4-way geometry, although we 
emphasize that the 4-way data are not treated as line-of-sight 
data, but as the fully integrated data that they are. Never- 
theless, with the 2-way data outnumbering the 4-way data by 
almost 5 to 1 , there is a strong effect present on the solution, 


which we can not avoid with our current data set. It might be 
possible to mitigate this effect slightly by re-weighting the 
data, but we did not pursue this, since we wanted to weight 
the data as close as possible to the weights used for 
SGMIOOh. In our benchmark test in section 3 such an orbital 
track effect was not present, but this is likely because there 


Angle Rstar-Main/radial 






Angle [deg] 


Figure 5. Geometry of the 4-way data over the SPA area, (left) The angle between the radial direction 
(taken as the ground track point of the main orbiter in lunar fixed coordinates) and the vector from the 
main satellite to the relay satellite are plotted, as well as (right) a histogram of counts of angles in 5° bins. 
Both plots show there is a substantial portion of the data with large angles, meaning less radial sensitivity. 
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the angle of line-of-sight with respect to the area of interest 
is different, and shorter arcs were used. 

[ 29 ] Finally, there is also more variation in the free-air 
adjustments in the southern part of SPA. This is because 
there are more data available there, due to the inclusion of 
2-way data around the south pole. These variations correlate 
with features such as Planck (57.9°S, 136. 8°E, 314 km), 
Schrodinger (75°S, 132.4°E, 312 km) and Poincare (56.7°S, 
1 63. 6°E, 319 km). 

[ 30 ] When comparing the hill solution (adjustment plus a 
priori field) to the a priori field only, the differences are more 
difficult to spot due to the different scales involved. When 
the two are displayed after each other in a loop (see auxiliary 
material), it becomes clear that the local adjustments modify 
the global a priori solution from SGMIOOh in such a way 
that several anomalies become less pronounced (e.g., the 
anomaly on the southern part of the ring around Apollo). 1 
The modifications of the features mentioned in the previous 
paragraph can be seen as well. They are more clear when 
viewed in this way, since their full signal is visible, rather 
than the adjustments only. The interior area of SPA is also 
slightly modified. 

[ 31 ] The adjustments themselves only show the differ- 
ences with the a priori field, and obviously from those alone 
one can not tell whether there really is an improvement or 
not. We therefore evaluated the local solution by looking at 
the 4-way data fit. And in the next subsections, we evaluate 
the local solution in a geophysical sense by comparing 
gravity and topography. 

[ 32 ] Data fit to the 4 -way data for the new solution is 
difficult to evaluate, due to the local nature of the adjust- 
ments: the concentration of the Slepian functions means that 
they should be evaluated only within the area, whereas the 
arcs used in our dynamical approach cover the whole Moon. 
In order to have a direct comparison, instead of using the 
initial arcs used for the processing, we used the short- 
arc approach as applied in the benchmark test over Mare 
Serenitatis, in this case on a few selected tracks of 4-way 
data over the SPA area. We also included one pass of 4-way 
data obtained on January 30 of 2009, that filled in a 
remaining gap over the farside [Matsumoto et al ., 2010]. 
During this pass, the main orbiter flew at an average altitude 
of 60 km instead of 100 km. We found that in general the 
local solution results in a slightly better 4-way data fit than 
the a priori model. As can be expected from the adjustments, 
the improvements are small: from 0.2 mm/s to 0.17 mm/s, in 
an RMS sense. For the 4-way data taken on January 30, 
2009, the improvement is from 0.24 mm/s to 0.18 mm/s. 
These values are smaller than the 4-way data fit values 
reported before by Matsumoto et al [2010], because they are 
for short arcs. Overall, the local solutions thus improve the 
data fit, albeit only slightly. The data fit improvements 
however are thought to be indicative yet inconclusive, since 
the tackiness in the solution might be an artifact of 
remaining orbital signals, causing the fit to be better while 
inducing spurious signal into the gravity solution. Hence we 
look at different ways of validation. 


Auxiliary material data sets are available at ftp://ftp.agu.org/apend/je/ 
201 lje00383 1 . Other auxiliary material files are in the HTML. doi:10.1029/ 
2011JE003831. 


4.3. Comparison With Lunar Topography 

[ 33 ] For the lack of direct and independent measurements 
of gravity on the planets, which could be used for evaluation 
and calibration of the models, the correlation between 
gravity and topography is often used to evaluate a gravity 
field model, with the idea that especially at the higher 
degrees of the spectrum uncompensated topography con- 
tributes directly to the measured gravity field, and thus 
higher correlations indicate an improvement. Here, we use 
the methods of Wieczorek and Simons [2005] to compute 
correlations and admittance between gravity and topography 
over the SPA area, and we also look at spatial correlations 
between gravity adjustments and craters on the lunar surface. 

[ 34 ] Figure 6 shows the correlation and admittance 
between various gravity models, and a topography model 
derived from SELENE LALT data [Araki et al , 2009]. The 
results also include a gravity model named SGM150, which 
is a preliminary 150th degree and order spherical harmonics 
expansion of the lunar gravity field, based on the SELENE 
data that were used for SGMIOOh, and historical tracking 
data coming from Lunar Orbiters 1-3, the Apollo 15 and 
16 sub-satellites (all satellites with non-polar inclinations), 
and all Lunar Prospector tracking data, including the 
extended mission data (which so far has not been included in 
the 100th degree and order SELENE models). LP extended 
mission data were processed in arc lengths of 2 days, using a 
data weight of 3 mm/s for Doppler data from January 1999 
until May 1999, and a data weight of 6 mm/s for Doppler 
data from May 1999 until the end of the mission. This 
division was chosen following Konopliv et al [2001], where 
different data weights were applied according to the location 
of the perilune of LP (over the nearside for the first half of 
the extended mission, and over the farside for the second 
half). SGM150 was furthermore estimated using a Kaula 
rule of 2.5 • 10 -4 // 2 , slightly stricter than what we applied 
to our local solution, but the same as other global solutions 
of the same resolution [Konopliv et al , 2001; Mazarico 
et al , 2010]. 

[ 35 ] The results in Figure 6 show that the local solution 
has the highest correlations with topography. The differ- 
ences are again not large, but they are pronounced for 
especially the degrees spanning / = 60-75. The tackiness in 
the adjustments (see Figure 4) does thus not affect the 
solution up to these degrees. After about degree 75, all 
solutions show a distinct drop in their correlation with 
topography, and this is likely because at this point the Kaula 
smoothing starts to become noticeable [Matsumoto et al , 
2010]. We used a spherical cap of 30° for the localization, 
in order to be safe from boundary effects in the local solu- 
tion. While 30° seems conservative, the correlations are 
always better for the local solution for larger cap sizes, but 
the differences between the gravity solutions become less 
pronounced. This means that, as often in local inversions, 
one should be careful with the solution at the boundaries. 
We also discuss this further in this section, when correlating 
topography and adjustments spatially. 

[ 36 ] At the higher end of the spectrum, beyond degree 90, 
SGM150 has better correlations than the local solution. This 
is probably due to the inclusion of the LP extended mission 
data into SGM150, showing that these data have an effect 
at the higher degrees (likely especially over the south pole, 
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Figure 6. (top) Localized correlations and (bottom) admittance over the SPA area, for a spherical cap of 
30°, for various gravity field models. 


as the 2-way data accumulate there), although for both 
models, correlations are dropping for the higher degrees. We 
discuss the LP extended data further in section 5. 

[ 37 ] The drop in correlations indicate that the solutions are 
losing their strength. The SELENE farside data assure that 
spherical harmonics models can be estimated up to degree 
and order 70 without the use of applying an a priori con- 
straint [Matsumoto et al ., 2010]. It is also beyond this degree 
that the correlations with topography start to drop. This 
might lead one to wonder whether the data allow a 150th 
degree and order solution. However, a 100th degree and 
order local solution does not show the improved correlations 
as the 150th degree and order solution shown in Figure 6 
does. This leads us to reason that even though the full res- 
olution of 150 by 150 might not be exploited, increasing 
the resolution at least assures an internal distribution of 
the signal over the model coefficients such that overall the 
solution improves up to about degree 75. This was also 
corroborated by a 120 degree and order local solution 
that we tested, that showed the same correlations as the 150 
degree and order local solution. 

[ 38 ] The admittance values (the transfer function between 
gravity and topography) are also shown in Figure 6. They 
show slightly different values in the same range of degrees 
as the difference in correlations did. All models are 
quite steady and level off at a value slightly higher than 


1 00 mGal/km. The local solution provides the most leveled 
and steady admittance value, at 104 mGal/km, whereas 
SGMIOOh reaches 102 mGal/km and SGM150 reaches 
99 mGal/km. We also point out that Figure 6 shows a second 
plateau for the admittance values for the local solution, 
between degrees 70 and 90, with a value of 85 mGal/km, 
which is not as pronounced in the other solutions. For 
comparison, we also included the admittance between 
LP150Q and the LALT topography, indicating how admit- 
tance estimates have changed with the inclusion of SELENE 
farside data. Admittance values can be used to estimate the 
lithospheric thickness [e.g., McKenzie and Nimmo , 1997; 
Crosby and McKenzie , 2005]. 

[ 39 ] To assess the correlations with topography spatially, 
we investigated the correspondence between the local grav- 
ity adjustments and craters on the lunar surface. We did this 
by plotting the crater database of Head et al [2010] (who 
produced a catalogue of impact craters >20 km in diameter 
on the lunar surface) on top of the adjustments evaluated 
between degrees / = 40-80, the result of which is shown in 
Figure 7. We chose this frequency band because between 
those degrees, the correlation with topography increased, as 
explained before. There is sometimes an offset between 
anomalies and craters, especially in the south-western part of 
the solution, which could be due to issues of data coverage 
or orbital errors, but for the larger part anomalies correlate 
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Figure 7. (left) Adjustments expressed as free-air gravity anomalies, in the frequency band that showed 
increased correlations with topography (/ = 40-80), and (right) differences between two local solutions 
using either data in a cap of 45° or 40°, to assess boundary effects. The 40° boundary is indicated 
by the dotted red line. Craters with sizes > 60 km, taken from Head et al. [2010], are indicated by 
black circles. 


with craters. For example, the outstanding adjustment at the 
rim of Apollo follows a trail of small-scale craters rather 
well, and several smaller craters show negative adjustments. 
Figure 7 thus shows that the adjustments correspond to 
surface features, as was also inferred from increased corre- 
lations with topography. 

[ 40 ] To assess the sensitivity of the adjustments with 
respect to the aforementioned boundary effects, we investi- 
gated the differences between two local solutions, based on 
localized functions in a cap of either 45° or 40° (and using 
data in those areas). These are also shown in Figure 7, again 
in the frequency band / = 40-80. As was the case with the 
results shown in section 3 when boundary effects were 
discussed, both solutions are affected by these effects. In 
addition, the difference plot in Figure 7 shows large differ- 
ences outside the 40° area because the 40° solution is 
undefined there. Differences between the solutions are of 
a smaller amplitude than the adjustments themselves, but 
nevertheless the interior of the solution area is affected, 
with the biggest differences within several degrees of the 40° 
boundary. Some of the differences shown in Figure 7 cor- 
relate with small craters, although those are mostly on or 
outside of the 40° boundary. Overall however, the interior of 
the SPA area seems to be stable between both solutions, 
which makes us reasonably confident of the adjustments. 

4.4. Bouguer Anomalies 

[ 41 ] One of the advantages of the Slepian functions is their 
direct interchangeability with spherical harmonics, which 
means that it is relatively straightforward to derive properties 
other than free-air anomalies from the estimated coefficients. 
Using the local solution and LALT topography, we com- 
puted Bouguer anomalies for the SPA area, assuming a 
crustal density of 2800 kg/m 3 and leaving out corrections for 


mare basalt. Bouguer anomalies for the local solution and 
SGMIOOh are shown in Figure 8. 

[ 42 ] The differences between the two solutions are more 
clear here than from the free-air anomalies, since the 
Bouguer anomalies have a wider range. Improvements in 
resolution for the local solution (while keeping in mind that 
the local solution is a larger expansion, in spherical harmo- 
nics) can be seen as certain features show more detail: for 
example, in the interior of SPA, around Apollo, Schrodinger, 
Planck, and in the vicinity of the Zeeman crater (75.2°S, 
133. 6°W, 190 km). Figure 9 shows the differences in 
Bouguer anomalies, with craters again plotted on top of 
the anomalies, and the aforementioned features are indicated 
as well. The differences are again most pronounced around 
the boundary areas, and while those are also influenced by 
boundary effects on the local solution, it should also be kept 
in mind that there is rougher topography at the edges of 
SPA (see Figure 1). Other differences correlate well with 
topographic surface features, as indicated by the agreement 
of anomaly differences with crater locations (although one 
should keep in mind that Bouguer anomalies were computed 
using topography data as input), especially in the southern 
part of SPA. This is possibly due to better data coverage 
there, despite trackiness induced by the 2-way data that has 
disappeared in the Bouguer adjustments. 

[ 43 ] From Bouguer anomalies, the compensation state can 
be derived. As expected, our results here do not change the 
overall compensation state of SPA as found by previous 
analysis [e.g., Zuber et al ., 1994], since our adjustments are 
relatively small with respect to the overall signal. While the 
extremes in Bouguer anomalies in Figure 9 vary between 
±200 mGal, the overall change is more modest, with an 
increase in the RMS of the Bouguer anomalies in the 
SPA area up to 289 mGal, from 283 mGal for SGMIOOh. 
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Figure 8. Bouguer anomalies over the SPA area from both the local solution and SGMIOOh. 


Understanding SPA’s compensation state may help in 
understanding how SPA contributed to the lunar orientation 
[ Garrick-Bethell and Zuber , 2009], and local differences 
within SPA itself may point to modifications to the crustal 
structure that occurred after formation of SPA [e.g., 
Neumann et al ., 1996]. In that respect, it is interesting to note 
that from Figure 8 a modification in Bouguer anomalies 
corresponding to the Zeeman crater (indicated in Figure 9) 
stands out. Two small circular negative Bouguer anomalies 
close to the south pole can be seen clearly in the local 
solution, while they were less visible (and not negative) in 
the SGMIOOh solution. The lower of the two corresponds 
to the Zeeman crater, while the upper of the two is at least 
partly covered by it. We point out that Zeeman is one of the 



Figure 9. Difference in Bouguer anomalies over the SPA 
area between the local solution and SGMIOOh, with craters 
with sizes > 60 km indicated. 


features where Yamamoto et al [2010] found exposures of 
olivine, which could lead to a different Bouguer signal, and a 
signature in the Moho uplift (see below). 

[ 44 ] Finally, we tested how our new solution modifies the 
estimate of the Moho topography beneath the SPA basin. 
We compare our results with those presented by Ishihara 
et al [2009] (which were based on a predecessor model of 
SGMIOOh, called SGMIOOg), and since those results were 
based on a global expansion of the gravity field, we trans- 
formed our local solution to an equivalent spherical harmo- 
nics expansion valid over the entire globe. We do this by 
simply generating free-air anomalies from the local solution 
within the SPA basin (in this case limited to a spherical 
cap of 32° radius to prevent boundary effects leaking into 
the solution), and from SGMIOOh outside the basin. We 
then estimate a global 100th degree and order model from 
these anomalies. This model has the same correlations with 
topography as the local model. There are some anomaly 
differences within the SPA basin when compared to the local 
solution, and this model does not take full advantage of the 
increased resolution of the local solution, but on the other 
hand, Ishihara et al [2009] used SGMIOOg up to degree and 
order 70, and we also find that up to this degree our local 
solution has increased correlations with topography. In other 
words, the loss of resolution should have only a small effect. 
Local crustal thickness analysis might be applied in order 
to make more use of increased resolution. 

[ 45 ] For determining the Moho topography, we expanded 
the models up to degree 80, using the same gravity inversion 
method as Ishihara et al [2009]. Figure 10 shows the Moho 
topography for different gravity field models beneath SPA. 
Their differences are also included, and from these, the 
modification on the southern rim of the Apollo basin that 
was also found in the free-air anomalies clearly stands out. 
It should be noted however that the feature disappears in the 
local solution, making the Moho structure south of the rim of 
Apollo a little smoother. The modification in the Zeeman 
crater that was found in the Bouguer anomalies can also be 
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Figure 10. Moho topography beneath SPA using different gravity field models and their differences. 


seen. From the differences only, it shows more uplift in the 
local model than it does for the global models, especially 
SGM150. However, our local solution does not seem to 
indicate a mantle plug associated with this area. Finally, the 
southern part of SPA is again affected most because there 
were more data there. 

5. Discussion 

[ 46 ] Here, we discuss the performance of the recovery 
method, with respect to specific choices we made, such as 
the choice of estimated parameters and the a priori model. 
We also briefly discuss the use of LP extended mission data. 

5.1. State Vector Effects 

[ 47 ] As mentioned in section 4.1, we formed the normal 
matrices per arc from the data over the SPA area, while 
including arc parameters (such as the adjustments to the 
initial states of the satellites involved, measurement biases 
and arc-dependent solar radiation pressure force model 
parameters). The integration of the orbit and computation of 
the measurements according the used models is then iterated, 
with the arc parameters adjusted per iteration (keeping 
the global parameters fixed) in order to fit to the tracking 
data, until convergence is reached. These arc parameters 
were then eliminated from the formed normal matrix by 
partitioning, during the aggregation of the separate normal 


matrices. We tested whether or not these parameters need to 
be included by also forming normal matrices using only the 
partial derivatives with respect to the gravity parameters. 
This is different from eliminating the arc parameters all 
together, since by choosing only the gravitational para- 
meters, the correlations between arc and common para- 
meters are not accounted for. We found however that the 
solutions using only the partials for the gravity coefficients 
did not show increased correlations with topography. Cor- 
relations in the range / = 60-75 dropped below those for 
SGMIOOh, and also those in the higher degrees (beyond 
100) were negatively affected. Gravity anomaly differences 
between these two local solutions (with and without the arc 
parameters included) are mostly confined to the southern 
part of SPA. The solution using only the gravity partials 
also shows more trackiness in its free-air anomalies. This 
indicates that there are still state vector effects present that 
need to be accounted for in the solution. Again, this likely 
affects the 2-way Doppler data most, since these data have 
coverage limited to the southern part of SPA. 

[ 48 ] During the partitioning step, the arc parameters are 
eliminated through back-substitution which involves the 
inverse of the sub-matrix containing the arc parameters [e.g., 
Kaula , 1966]. We use all the data in the two-day arcs to 
converge them, but using only the data over the area of 
interest to form the normal matrix might lead to numerical 
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instabilities in the inversion of especially the arc parameters 
part when the normal matrices are aggregated. This is done 
with the SOLVE software, which allows for a dynamical 
suppression in the Cholesky inversion of parameters that 
turn out to have poor observability (indicated by very small 
or even negative numbers on the diagonal). Inspection of the 
results has shown that this is sometimes the case, mostly for 
certain bias-parameters that were included, but occasionally 
also for state vector parameters. Since we wanted to include 
all available 4-way data, we did not exclude those normal 
matrices where this was the case from the solution. 

5.2. Influence of the a Priori Gravity Field Model 

[ 49 ] In order to test the sensitivity of the method with 
respect to the a priori model, we processed the same data 
using LP100K [Konopliv and Yuan , 1999] as the a priori 
gravity field model. In this case the a priori model is inde- 
pendent of the SELENE data, as LP100K is a 100th degree 
and order spherical harmonics model based on LP tracking 
data (and as such, a precursor to the LP165P and LP150Q 
models). This test is similar in extent to the benchmark 
test presented in section 3, in that we process the data with 
an independent model to test whether the information can 
be extracted from the data. For LP100K, the gravity field 
over the SPA region is much smoother than that described 
by SGMIOOh, and a basin such as Apollo is not circular 
in shape. Correlations with topography over the SPA area 
are also much lower than those for the SGMIOOh and 
SGM150 models, with a distinct low around degree 60 
where the correlations are close to zero. The SELENE 4-way 
data show large residuals with respect to LP100K [Namiki 
et al ., 2009]. 

[ 50 ] We found that our local solution based on LP100K as 
a priori model resolves the structure within the SPA basin 
rather well. The Apollo, Planck and Ingenii features stand 
out much clearer in the local solution than they do in 
anomalies computed from LP100K. Yet there is also some 
orbital trackiness visible in the solution, and while correla- 
tions with topography improve drastically over those for 
LP100K, up to correlations of 0.89 for the range / = 40-70, 
they do not reach the levels as obtained by the global or local 
SELENE models, which is about 0.93 for the same range. 
This means that we can only adjust the a priori model to 
some extent, and that, unsurprisingly, we get the best solu- 
tions if we start with a model that already describes the data 
well. One should also keep in mind that the dynamical 
approach in processing the tracking data here assumes line- 
arity, and that global spherical harmonics models are often 
determined from multiple iterations over the data. Due to the 
local nature of our representation of the gravity field, which 
should only be evaluated within the area of interest, we did 
not iterate our solutions, since that would require processing 
arcs which also cover the rest of the Moon, outside of SPA. 

5.3. LP Extended Mission Data 

[ 51 ] LP extended mission data are especially interesting 
for local nearside analysis, and they have been the focus of 
previous local and regional analyses. And since we have also 
included SELENE 2-way Doppler data beyond the south 
pole, LP extended mission data might be included as well. 
We processed the LP extended mission data in the same 
way as they were processed for the SGM150 model (see 


section 4.3). We again selected only data over the SPA area 
when adding them to our local solution. Since these data 
are not included in SGMIOOh, they showed large residuals, 
and this has a profound effect on the solution: it shows large 
unrealistic gravity anomaly fluctuations. 

[ 52 ] When we compared various gravity models in terms 
of their correlations with topography (see section 4.3), we 
noted that the SGM150 model showed increased correlations 
at the higher degrees, likely coming from the influence of LP 
extended mission data. In order to bypass the large residuals 
in the LP extended mission data when we process them with 
SGMIOOh, we recomputed the local solution starting from 
SGM150 instead (which included the LP extended mission 
data), using both SELENE and LP extended mission data. 
The SELENE-only solution shows similar free-air anom- 
alies and correlations as the local solution starting from 
SGMIOOh. When we include the LP extended mission data 
however, correlations with topography unfortunately drop 
again, and orbital trackiness starts to appear in the solutions. 
If we down-weight the LP extended mission data when 
combining the normal matrix with that of the SELENE data, 
the solution improves, but it does not lead to for example 
increased correlations when compared to the SELENE data- 
only solution. Sensitivity issues as discussed in section 4.2 
might especially play an important role here, and thus the 
LP extended data were not included into our solutions. 

6. Conclusions 

[ 53 ] We have estimated adjustments to a global lunar 
gravity field model over a limited area, by representing the 
gravity field by means of localized basis functions. We 
applied Slepian basis functions, which are represented 
by linear combinations of standard spherical harmonics. 
This makes them relatively straightforward to use in com- 
bination with existing data processing strategies based on 
spherical harmonics. We process tracking data in the stan- 
dard dynamical way, which means that we follow an inte- 
grated approach where the local adjustments are estimated 
from the tracking data residuals directly (as opposed to 
numerically differentiating the Doppler residuals to obtain 
accelerations). This is especially useful for tracking data that 
consist of the accumulated Doppler shift between different 
links, such as the SELENE 4-way data. We generate the 
design matrix (which describes the linearized relationship 
between the residuals and the estimated parameters) using 
only data over the area of interest. From these reduced 
design matrices we create normal matrices and aggregate 
them in the standard way. The normal matrix is then trans- 
formed from spherical harmonics into Slepian functions, and 
we only estimate those function coefficients with a concen- 
tration factor larger than a pre-set value, which means that 
only functions that have power in the area of interest are 
estimated. The estimation is done in a least-squares sense. 

[ 54 ] A benchmark test for the gravity recovery method 
was conducted, using Lunar Prospector (LP) extended 
mission tracking data over the Mare Serenitatis area on 
the nearside of the Moon. The data were processed with a 
pre-LP lunar gravity field model (GLGM-2) that shows 
smoothed gravity anomalies over Serenitatis when compared 
with LP-based lunar gravity field models. The test showed 
that the method as implemented can extract the high 
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resolution information from the tracking data, as we found 
improved correlations with anomalies from a high resolution 
global lunar gravity field model using LP data. 

[ 55 ] We then estimated adjustments to the 100th degree 
and order SGMIOOh model for the SPA area, using Slepian 
functions in a resolution equivalent to degree and order 150 
in spherical harmonics over an area with a spherical cap of 
45° radius, centered at (191. 1°E, 53.2°S), resulting in the 
need to estimate 4,446 Slepian function coefficients (as 
opposed to the 22,797 needed for a full global spherical 
harmonics model). The use of a limited set of data (only 
those over the area of interest) also means that we can create 
our local solutions much faster than a global solution of 
equivalent resolution. We applied a Kaula rule of 3.6 • 1 O -4 / 
/ 2 to this solution to prevent large fluctuations in gravity 
anomalies. We used SELENE tracking data, including 4- 
way data over the SPA area, and 2-way data slightly beyond 
the south pole of the Moon. Even though we estimate fewer 
coefficients, the resolution of the model and the satellite 
altitude of 100 km are such that smoothing is required. An 
investigation into the effects of choosing data only over the 
area of interest showed that the solutions are generally 
affected up to a few degrees into the area of interest. 

[ 56 ] The adjustments are small in general, between plus 
and minus 68 mGal with an RMS of 17 mGal. Foremost, this 
indicates that the a priori model SGMIOOh captures the 
information in the 4- way data well. Our new local solution 
shows increased correlations with high-resolution lunar 
topography, and admittance values that are slightly different 
and more leveled when compared to those from global 
models using the same SELENE data. Free-air anomaly 
adjustments, evaluated in the frequency band / = 40-80 
(where the increase in correlations with topography was 
concentrated) correlate with craters on the lunar surface. 
Since our local solution is easily transformed into spherical 
harmonics, we also computed Bouguer anomalies. The local 
solution shows an increased resolution inside the SPA area, 
and differences with Bouguer anomalies from the global 
models correlate well with topographic surface features. The 
overall compensation state of SPA as found by previous 
analysis is not changed by our results, since we mostly 
adjust the global lunar gravity field models at their short- 
wavelength components. The Moho structure beneath the 
SPA basin is slightly modified by our solution, most notably 
at the southern rim of the Apollo basin, where a previously 
outstanding feature now appears smooth. We also found 
modifications in Moho structure and Bouguer anomalies for 
the area corresponding to the Zeeman crater (where olivine 
exposures were found). 

[ 57 ] We demonstrated that the method can complement 
global solutions by refining them over areas where data are 
more dense. This might be especially of interest for future 
planetary missions such as GRAIL [Zuber et al . , 2008], 
where the lunar gravity field will be mapped from low-low 
satellite tracking data. While it is outside of the scope of this 
work to do a full analysis, such a data type is highly sensitive 
to local gravity field variations, and thus suitable for a local 
method such as presented here. 

[ 58 ] The local solution over the SPA area that we pre- 
sented here is made publicly available as a set of normal- 
ized spherical harmonics coefficients through the RISE Data 


Archive Web site, which can be found at http://www.miz. 
nao . ac .jp/rise-pub/en. 
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